Assigned Reading: Geocomputation with R Chapter 2.

What is a Geographic Information System?

What is a GIS

• Geographic Information System is a system for producing, managing, displaying, and analyzing geographic information.
• Spatial Data Science is an emergent field which utilizes data science approaches in a GIS, and is a natural extension of a GIS

A Visualization of a GIS, by Anne Sexton

A Brief History

• Dr. John Snow suspected Cholera was spread by water
• 1854 outbreak of Cholera in the Soho neighborhood of London kills 616 persons
• Snow used both mapping and statistics to identify the water source and stop the outbreak.
• Essentially founded both Epidemiology and GIS, and stopped the outbreak in one stroke

The Broad Street Pump, A Cholera Outbreak in London, By: John Snow and digitized by National Geographic

https://www.theguardian.com/news/datablog/2013/mar/15/john-snow-cholera-map

• Since then GIS waited to flourish until computational ability increased dramatically in the 90’s

What is Spatial Data Science?

• The application of Geographic insights, and geospatial analyses to big data sets
• Including spatial terms in statistical models
• Bringing Geographic information to data science questions
• Developing your own spatial products and pipelines

Why Use R as a GIS and for Spatial Data Science ?

• Work Flow Automation
• Rich Ecosystem (packages, functions, code-sharing, etc.)
• Reproducible
• Self Documenting
• Computationally Efficient
• Parallel Processing/HPC interfaces

Geodesy

• ‘Geodesy is the science of accurately measuring and understanding the Earth’s geometric shape, orientation in space, and gravity field.’ - NOAA

Earth is not a perfect sphere

• Circumference at equator: 40,075 km (24,901 mi)
• Circumference along meridians: 40,009 km (24,860 mi)

Blue Marble 2012, by Suomi NPP

The Earth can be represented as an Ellipsoid

• A simple 3d geometric shape
• An ellipsoid is a slightly to greatly ovaliform shape
• Modeling the Earth as an ellipsoid increases the accuracy of point locations

Ellipsoids

However, the Earth’s surface is not smooth - Geoid

• Types of models to represent planet Earth
• Include gravity, excluding winds and tides
• Highly accurate since application of GPS technology

Left: Geoid Cross Section. Right: IGCM Geoid

By International Centre for Global Earth Models (ICGEM) - http://icgem.gfz-potsdam.de/vis3d/longtime / Ince, E. S., Barthelmes, F., Reißland, S., Elger, K., Förste, C., Flechtner, F., Schuh, H. (2019): ICGEM – 15 years of successful collection and distribution of global gravitational models, associated services and future plans. - Earth System Science Data, 11, pp. 647-674,DOI: http://doi.org/10.5194/essd-11-647-2019., CC BY 4.0, https://commons.wikimedia.org/w/index.php?curid=81462823

Geodetic Datums

• Reference frame established to represent locations within the frame.

• Historically these were locally focused and based on Geoids.

• A datum can serve either horizontal (X & Y) or vertical (Z) features.

• Components:

  • reference ellipsoid or geoid
  • origin point (from which measurements run)
  • control points very strictly measured from the origin,

Other locations are then measured in relation to the control points.

Number one City Datum. by: Cosmo1976

Control Points Number one City Datum

Geodetic Coordinatees

“The geodetic (or geographic) latitude is the angle between the equatorial plane and the normal (vertical) to the ellipsoid surface at the considered point.” - Proj

• phi = geodetic latitude (north/south)
• lambda =longitude (east/west)
• h = ellipsoidal height
• N = Normal (A plane at a right angle to the surface of the ellipsoid)

Angles on Ellipsoid and Geodetic Coordinates. by: Peter Mercator

By Peter Mercator - Own work, CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=17717979

Coordinate Notation

• Trigonometric
  • Degrees Minutes Seconds (DMS), e.g. 42°03’27.7” N
    • Sexagesimal/base 60 (think: minutes in an hour)
  • Decimal Degrees (DD), e.g. 42.05759 N
    • Decimal Fractions of a Degree (portion of 360/2)
• DMS; seldom used in digital formats
• DD; almost exclusively used

\[ \text{Decimal Degrees} = \text{Degrees } + \frac{\text{Minutes}}{60} + \frac{\text{Seconds}}{3600} \] \[ \text{Latitude of Tech} = 42°03'27.7" N \]

\[ 42 + (\frac{03'}{60}) + (\frac{27.7"}{3600}) \]

\[ 42 + 0.05 + 0.007694 = 42.05759 \text{ Decimal Degrees} \]

• Universal Transverse Mercator (UTM)
  • Divides the world into 60 zones
  • Flattens each zone
  • Measures distances in Meters

UTM - Common for field work, planar measurements of meters.

UTM Zones of the Conterminous USA. data by ESRI

Coordinate Reference System

• System for specifying a location on Earth’s surface
• Composed of:
  • a model of earths shape (Geoid or Ellipsoid)
  • geodetic datum
  • generally also a projection

Geographic Coordinate Reference System. by: Anna Krystalli

Problems with working on flat surfaces

3D to 2d does not work well. But we have worked around this. sorta

Orange Peel. by: Nathan Belz

Geographic Coordinate System

• Location on a three-dimensional model of Earths Surface

Projected Coordinate System

• Location on a two-dimensional model of Earths Surface

Major Map Projections

Many different ways to create two dimensional maps – thousands of projections. None are perfect. There are four main types of projections, each with their own strengths and examples. One of the most common examples of each of these is included in the table below.

• Equal Area: Lambert Cylindrical equal-area
  • Pro: No distortion of area near equator (here)
  • Con: Distorts area at the Polar regions
• Equal Distance: Equi-rectangular (plate carrée projection)
  • Pro: Looks good in mapping applications
  • Con: Distorts both shape and directions
• Conformal:
  • Pro: Boundaries are accurate
  • Con: Distorts Polar area
• Compromise:
  • Pro: Sensitive to both area and direction
  • Con: Sensitive to both area and direction

Map Projections. by: Daniel Strebe

Geodesy Takeaways

Unless you really focus on GIS, hardly anything I have said this lecture will matter to you.

• There are various models (geoids & ellipsoids) to represent the shape of the earth
• Different datums are used to represent different parts of the earth. This is in part due to legacy effects.
• You will almost always use WGS 84 (based on a ellipsoid – which is fit through a special model of earths gravitational fields a geoid) NAD83 (based on a ellipsoid) for a geographic coordinate system. These +/-1 m from each other across much of North America. More useful than a geoid.
• You will usually want to use a UTM grid based on WGS for or State Planes based on NAD83 for projections.
• Different coordinates notation systems are used, focus on Decimal degrees.

An Introduction to Geographic Data Models

• Vector Data Model represents discrete features on the planet using geometries such as: points, lines, and polygons.
• Raster Data Model represents (usually) continuous features on the plant using continuous surfaces, like a tile.

Vector data tends to only include features of interest, e.g. bodies of water; whereas a Raster will include, explicitly, the absence of features (e.g. both water and terrestrial areas).

Vector and Raster Data Models

R Data Types Reviewed

Remember that different data types take up different amounts of memory in pretty much all software systems.

When the values in raster cells are stored as integers as opposed to floats (with decimal points) they take up only 0.5 as many bytes.

When the values in raster cells are stored as integers as opposed to words they take up only 0.5 as many bytes.

The size of the our empty raster is 0.01416 MB

The size of the data for our raster as a character vector is 0.11156 MB

The size of the data for our raster as a numeric (integer !) vector is 0.055744 MB

The size of the data for our raster in matrix form is 0.055912 MB

The size of the data for our raster in dataframe form is 0.055912 MB

The size of our complete raster is 0.069856 MB

Vector Data in R

• Vector Data Model represents discrete features on the planet using geometries such as: points, lines, and polygons.

Simple Features - Standards

• Open Geospatial Consortium ISO 19125-1:2004: Currently adhered to in ESRI, used in GDAL.
• Features have geometries describing their locality on Earth, and properties described by attributes.
• If the geometry of a feature is a polygon, it is composed of points, connected by straight lines.
• Lines composing polygons cannot intersect

Simple Features 1 - Attributes

• Attributes of the feature, theoretically devoid of spatial context.
• Described in text, numbers, stored in a data frame type object.

Simple Features 2 - Coordinates form a Point

• all geometries are composed of points.

• points only require two coordinates, X & Y.

• Y = Latitude (necessary), X = Longitude (necessary)

• Z = Elevation (somewhat uncommon)

• M = Time or Uncertainty of Measurement (rather uncommon)

Simple Features 3 - Points are/form a SF Geometry (sfg’s)

• SFG is the spatial topology associated with a feature

• POINT - A true dimensionless 1 dimensional point

• LINESTRING - Sequence of points connected by strings.

• POLYGON - Sequence of points connected by strings, which close upon themselves.

• i.e. the first and last point of the polygon are the same.

Simple Features 4 - A Geometry can be combined to form Geometries (sfg’s)

• Two or more Points, and sets of Linestring, and polygons can be the geometries of a single feature
• Multi(POINT), multi(LINESTRING), multi(POLYGON) - Forming collections of geometries
• GEOMETRYCOLLECTION - A mixed set of these three or more geometries in the same geometry

Simple Features 5 - Geometries and Spatial Information form an SF Collection (sfc)

• A list column of S3 type containing certain parameters regarding the Geometry/Geometries

• Coordinate Reference System (‘crs’)

• Precision (‘precision’)

• Bounding box (‘bbox’)

• Number Empty (‘n_empty’)

Wherein precision refers to the precision of the coordinates forming the geometries. the bbox is a rectangle which includes all pairs of coordinates in the geometry. N_empty, denotes whether the geometry is missing coordinates in a position.

Simple Features 6 - Recapped.

Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -87.67607 ymin: 42.05663 xmax: -87.67607 ymax: 42.05663
Geodetic CRS:  WGS 84
Precision:     2 

sfc_POINT of length 1; first list element:  'XY' num [1:2] -87.7 42.1

• We now finally have both attributes and coordinates which are forming a point geometry.
• This is a Simple Feature.

sp

• Predecessor to the ‘SF’ R package
• Many spatial statistics programs still only take sp objects as input
• A good introduction to S4 objects; popular in new spatial packages too!

sp - History

• Released in 2005

• First package to hold all major vector types of geometries

• Unified Spatial class allowed for support in many spatial statistics packages

• Improved mapping of spatial objects

• Formed the centerpiece of spatial statistics in R for over a dozen years.

sp - Spatial Classes for Topology

• SpatialPoints
• SpatialLines
• SpatialPolygons
• SpatialGrids

Similar in theory to the simple feature geometries we learned about…

sp - Structures of The Spatial* Class - Topology Only

• Bounding box
• Coordinate Reference System
• …Coordinates!

sp - Attributes - DataFrame

• Attachment of attributes, in a data frame, to a Spatial topology can create:

• SpatialPointsDataFrame

• SpatialLinesDataFrame

• SpatialPolgonsDataFrame

• SpatialGridsDataFrame

• Each S*DF is accessed the same way e.g. SPDF@data[[‘col_name’]]

We see here that SF contains all of the spatial information in a Simple Feature Collection which can be tidily tucked away in the geometry list column, in SP the analogue to the SFC is the object itself.

SP is mostly geometry with a data slot, SF is mostly data frame with a geometry column.

sp Overview of a Spatial*DataFrame

• Data frame is held in a different slot from the geometry and topology
• Data frame columns may be subset
• Data frame columns accessed via object@data[[‘colname’]] indexing

sp - Recapped.

• If you need to use spatial statistics, you will come across these objects.
• Don’t sweat them too much; you can always convert from and to sf to run certain analyses

Raster data in R

• Reminder: the Raster data Model is a way to represent (continuous) features on the planet using grids
• Excellent format for storing and sharing data
• Divides space into continuous grids
• Raster data sets are often distributed in tiles

Western Plant Predictors. Note each cell is the extent of a single raster tile, each tile contains cells.

Raster Components

• Bounding Coordinate(s)
• Cell Resolution (Size of cell)
• Dimensions (No of cells in rows and columns)
• Coordinate Reference System (CRS)

Note that both cell sizes, and the resolution of the values in cells can, within reason, be converted to finer and coarser resolution. We will discuss these types of calculations next class.

Example Raster 1 Create Frame

Rasters tend to confuse people. We are going to create our own example here before we talk about them much more.

Fairy shrimp live in bodies of water which dry out in late Spring, and refill in early Fall. We seek to determine where suitable habitat for fairy shrimp are in the Great Basin. In order to do so, we will use satellite imagery to detect valleys which fill with water in Fall, and dry out in late Spring.

Suprise Valley by: John Glen

empty_raster <- raster(
  
  # rasters have 4 bounding edges
  # Here we define each 'corner' of the raster'
    xmn = 697129.7, 
    xmx = 811775.7,
    ymn = 4388466,
    ymx = 4502382,
    
  # Here we set the number of cells 118*118
    nrows = 118, 
    ncols = 118,
    
    # we do NOT manually specify the cell size here.
    
    crs = "+proj=utm +zone=10 +datum=WGS84", 
    # set the rasters Coordinate Reference System
)

In the code above, we are going to define all four of the essential components of a raster. Clearly, we are defining three explicitly (Bounding Coordinates, Dimensions, CRS), and the remaining one implicitly (Cell resolution).

We can imagine that the raster we are currently creating looks like this. A frame without content. We can see where the bounding coordinates are,

Example Raster 2 Set Values

raster_matrix <- matrix(rast_vals_num, # # fill matrix with values, 
                       nrow = empty_raster@nrows, # create matrix same dimensions.
                       ncol = empty_raster@ncols, # create matrix same dimensions
#ensure filling of matrix goes from upper left to lower right.
                       byrow = T) 

raster_matrix[90:118,112:118] <- 1 # fix the clipping image at edge. 

I think it is easiest to imagine that the values within a raster are in the shape of a matrix. Hence, using the code above we could take a vector of values, and bend them, so that each position within the vector matches up with the beginning of a new row.

example_raster_dec <- setValues(empty_raster, raster_matrix)

Example Raster 3

The width of each raster cell is: 971.57627 meters
The height of each raster cell is: 965.38983 meters

This example_raster_dec contains: 13924 elements

Raster Package - in R!

• Does not need to load all files into active memory at once
• Many functions based on functions in ‘base’ R.

Attributes

• Cells
• Values

Raster Layer

A single raster file.

Raster Brick

As I have mentioned rasters are often generated from satellite imagery; we will discuss rasters which are not developed this way next lecture. Historically most pictures are imaged via the use of three bands. These having spectral values of Red, Green, and Blue (RGB) associated with them. For example our .tif file, is composed of three bands. Note that a Rasterbrick is most often used for loading in these types of imaging data, which can then be processed to form a more typical ‘raster’ dataset.

Rasters are often built from satellite imagery. Historically most pictures are split into three copies, each one having spectral values of Red, Green, and Blue (RGB) associated with them. For example our .tif file, is composed of three layers. Note that a Rasterbrick is most often used for loading in these types of data.

Bricks are very important for processing and classifying raw image data. While we have three bands here, nowadays LIDAR equipped with Hyperspectral sensors are likely to have many more bands up to a couple hundred depending on the application.

One important technical point to note with the Rasterbrick is that each band of the brick loads from the same individual file, for example picture. And that bands are not typically combined from different picture sources.

Raster bands also do not need to be held in memor, allowing one to work through large amounts of them.

Raster Stack

In general these are rasters which have been classified and we want to extract values from or run calculations with. Now what is great about layers, is that we can do one big thing bricks cannot do, we can load in many layers from many files and create a stack of attributes we are interested in studying on the fly.

The main use of Raster Stacks is to hold layers of similar themes.

Two common examples: 1) Each raster layer in the stack is a variable from a different month, e.g. mean monthly temperature from January -> December (12 layers per stack) 2) Each raster layer in the stack is a variable of interest in an analysis, e.g. yearly mean temperature, mean precipitation, etc. from which we want to extract values.

Does not need be held in memory

Terra

• Developed by the same team as the Raster Package
• Functionally virtually identical to Raster, but calculations run more quickly ! :-)
• Rasterbricks/layers no longer need specification, all objects are Spatstats ?
• Will supercede Raster, but see points 1 & 2.

Assignments for the Duration of the Spatial Data Science Module:

For next Lab:

Please install these packages (if you have not done so already):

install.packages("sf", "raster", "terra", "sp", "tmap", "leaflet", "ggmap")
# optional packages for using parallel processing at a step
# (in our example it won't actually speed up the analyses 
# they actually may slow them down!)
install.packages('snow','parallel')

Download the labs .R script from the course website.

If you are interested in how Drone/Lidar data are collected please check out the ‘RMBL Spatial Data Science Webinar Series’:

• https://github.com/ikb-rmbl/SpatialDataScienceWebinars2020
• Collecting UAS Data (Video): https://youtu.be/Pq8btEZRCvM (1.25 hours)

For next Lecture: Assigned Reading: Chapter 3 of Spatial Data Science

Future Bonus SDS Office Hours Wednesday Night at 5:00 - 6:00.

Notes on this Lab My lecture notes are in the R script I used to generate all of the novel figures for this presentation. Likewise this presentation is an .HTML file and can be launched from your computer (it was rendered directly from R using the script).

Works Cited

Krystalli, A. https://annakrystalli.me/intro-r-gis/gis.html Accessed 01.20.2022

https://geocompr.robinlovelace.net/spatial-class.html Accessed 01.09.2022

Hijman, R. 05.12.2019 ‘The raster Package’

https://rspatial.org/raster/RasterPackage.pdf Accessed 01.09.2022

https://www.neonscience.org/resources/learning-hub/tutorials/dc-multiband-rasters-r Accessed 01.19.2022

Pebesma, E. https://r-spatial.github.io/sf/articles/sf1.html Accessed 01.10.2022

https://proj.org/operations/conversions/geoc.html Accessed 01.20.2022.

https://rspatial.org/raster/spatial/8-rastermanip.html Accessed 01.11.2022

https://en.wikipedia.org/wiki/Open_Source_Geospatial_Foundation Accessed 01.09.2022

NOAA. What is geodesy? National Ocean Service website, https://oceanservice.noaa.gov/facts/geodesy.html, 1/25/17.

Geocomputation with R. Lovelace, R., Nowosad, J., Muenchow, J. 2022-01-25.